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Abstract 

Structural reliability and decomposition techniques have recently proved to 
be appropriate tools for solving robust uncertain mixed-integer linear pro¬ 
grams using ellipsoidal uncertainty sets. In fact, its computational perfor¬ 
mance makes this type of problem to be an alternative method in terms of 
tractability with respect to robust problems based on cardinality constrained 
uncertainty sets. This paper extends the use of these techniques for solving 
an adaptive robust optimization (ARO) problem, i.e. the adaptive robust 
solution of the transmission network expansion planning for energy systems. 
The formulation of this type of problem materializes on a three-level mixed- 
integer optimization formulation, which based on structural reliability meth¬ 
ods, can be solved using an ad-hoc decomposition technique. The method 
allows the use of the correlation structure of the uncertain variables involved 
by means of their variance-covariance matrix, and besides, it provides a new 
interpretation of the robust problem based on quantile optimization. We also 
compare results with respect to robust optimization methods that consider 
cardinality constrained uncertainty sets. Numerical results on an illustra¬ 
tive example, the IEEE-24 and IEEE 118-bus test systems demonstrate that 
the algorithm is comparable in terms of computational performance with 
respect to existing robust methods with the additional advantage that the 
correlation structure of the uncertain variables involved can be incorporated 
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straightforwardly. 


1. Introduction 

Robust optimization techniques were developed to deal with the prob¬ 
lem of designing solutions that are immune to data uncertainty by solving 
equivalent deterministic problems, i.e. robust counterparts (RC). The hrst 
published work about these techniques is proposed by Soyters [1], which 
allowed uncertain parameters to take their worst possible values within a 
given interval. However, although using this method the RC of any linear 
mathematical programming problem remains linear, Soyters’s solutions are 
excessively conservative, and [2, 3, 4] proposes the use of alternative uncer¬ 
tainty sets, the ellipsoidal uncertainty sets based on the Mahalanobis dis¬ 
tance. This alternative method allows controlling the grade of conservatism, 
and includes straightforwardly the correlation structure of the uncertain pa¬ 
rameters involved by means of the variance-covariance matrix. In contrast, 
the RC derived from any linear mathematical programming problem using 
ellipsoidal uncertainty sets results in a conic quadratic problem. This feature, 
which is not problematic for those cases where interior point algorithms are 
applicable because they have been proved to be solvable in polynomial time, 
is not particularly attractive for solving robust linear discrete optimization 
models, as pointed out by [5]. It is precisely [5] the one that proposes the 
use of an alternative uncertainty set that allows the RC problem to remain 
as linear, even if binary variables are involved. Instead of allowing all uncer¬ 
tain parameters to take their worst possible values within the given interval 
such as [1], it only allows a pre-established maximum number of parameters 
T to do so. This alternative formulation also provides a robust solution in 
terms of probability of infeasibility. Recently, [6] proposes an alternative 
and decomposable solution technique inspired on structural reliability First- 
Order-Second-Moment (FOSM) methods, which allows reaching the optimal 
solution of the RC problem using ellipsoidal uncertainty sets by solving two 
kind of problems within an iterative scheme: one master problem of the 
same size as the original RC problem, and one subproblem of considerable 
lower size for each uncertain constraint. The method proposed by [6] has the 
following advantages: 

1. It improves tractability with respect to interior point algorithms, spe¬ 
cially for large-scale problems and those including binary decision vari- 
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ables. 

2. In addition and due to its relationship with reliability-based structural 
techniques, it allows to calculate exact bounds for the probability of 
constraint violation in case the probability distributions of uncertain 
parameters are completely dehned by using hrst and second moments 
(mean and variance-covariance). 

3. In terms of tractability, it constitutes an alternative formulation, spe¬ 
cially for problems including binary variables, to that proposed by [5] 
using cardinality constrained uncertainty sets. 

Note that the formulation proposed by [6] was tested and proved to be 
successful on mixed-integer linear mathematical programming problems with 
hard constraints that must be satished for any possible realization of the un¬ 
certain data. This type of problems is relevant for optimal design, however, 
problems with recourse [7, 8, 9] where the decision-maker must take some 
decisions before discovering the actual value of the uncertain parameters and 
having the opportunity to take further action once the uncertain param¬ 
eters have been revealed, are more common on engineering problems and 
can not take advantage of the features of the method proposed by [6] in its 
present form. This framework is also known as Adaptive Robust Optimiza¬ 
tion (ARO). 

A clear example of the practical importance of robust linear optimiza¬ 
tion with recourse (ARO) is the Transmission Network Expansion Planning 
(TNEP) problem. It aims to resolve the issue of how to expand or reinforce 
an existing electricity transmission network so that adequately serves system 
loads over a given horizon. The main difficulty of this problem is to take 
decisions under the great amount of uncertainty associated with i) demands, 
ii) renewable generation, such as wind and solar power plants, and hi) equip¬ 
ment failure. There are several examples of the successful application of 
robust optimization in TNEP problems based on cardinality constrained un¬ 
certainty sets proposed by [5], such as, [10, 11, 12, 13, 14, 15]. These works 
use three-level formulations where hrst level minimizes the cost of expan¬ 
sion ([13] also minimizes the maximum regret), the decision variables for this 
level are those associated with construction or expansion of lines, ii) second 
level selects the realization of the uncertain parameters that maximizes the 
system’s operating costs within the uncertainty set, the variables related to 
this level are the uncertain parameters, i.e., generation capacity and demand, 
and iii) in the third level the system operator selects the optimal decision 
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variables to minimize operating costs for given values of first and second level 
variables. 

This paper proposes a double iterative solution approach to the three- 
level mixed-integer optimization problem associated with TNEP. The idea 
is to replace cardinality constrained by ellipsoidal uncertainty sets, which 
allows incorporating the correlation structure of the uncertain parameters 
straightforwardly by means of the variance-covariance matrix. Inner level 
problems (subproblems) are solved throughout an iterative method inspired 
from structural reliability and decomposition techniques [6]. Note that, due 
to the way the uncertainty sets are constructed, the subproblem does not 
contain binary variables. Meanwhile, first level problem (master problem) 
uses the constraint-and-column generation method (proposed by [16] ) ini¬ 
tially used in this type of problems by [13] and [14], and also utilized by 
|16|. 

The proposed method has the following advantages which might make it 
attractive for practical use: 

1. Computational performance is comparable with the method proposed 
by [15] based on cardinality constrained uncertainty sets. Note that 
this method proved to be the most efficient method to date for solving 
this type of problem. 

2. The subproblem can be solved using two nested linear optimization 
problems: i) the first one takes the optimal operating decisions for 
given realizations of the uncertain parameters involved, and ii) the sec¬ 
ond one updates the worst possible realization of the uncertain param¬ 
eters in terms of operating costs. This second problem, under certain 
circumstances, have analytical solutions and otherwise can be solved 
efficiently by using interior point algorithms. 

3. The approach proposed in this paper requires the definition of the pro¬ 
tection level, which establishes approximately the quantile of the uncer¬ 
tain objective function to be optimized. Besides, the proposed method 
would not require the definition of bounds for the uncertain parameters, 
which might be of interest for certain cases. 

The remainder of the paper is structured as follows. Section 2 describes 
the adaptive robust formulation of the TNEP problem in compact matrix 
form, and discusses the solution approach proposed by [15] inspired in the 
works by [12, 13, 14]. In Section 3 the definition of the uncertainty set and 
the description of the solution approach are presented. Numerical results 
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for different networks are given in Section 4 and compared with those ob¬ 
tained using the method proposed by [15]. Finally, the paper is concluded in 
Section 5. 


2. Transmission network expansion planning problem: ARO com¬ 
pact formnlation 

Note that for clarity in the exposition, we adopt the same nomenclature 
used by [15]. Thus, the robust TNEP problem can be written in compact 
matrix form as the following three-level mathematical programming problem: 


Minimize 

X 


cFx + 


Maximum 

deD 


Minimum y 
y E il{x, d) 


subject to 


cP^x < n 

X E {0,1}, 


( 1 ) 

( 2 ) 

( 3 ) 


where x is the vector of hrst stage binary variables representing the invest¬ 
ment vs no investment in reinforcing or building new lines, c is the investment 
cost vector, d is the vector of second stage continuous variables representing 
the random or uncertain parameters, i.e. generation capacities and demands, 
D is the uncertainty set, b is the vector including operating costs, and y 
is the vector of third stage continuous variables referring to the operating 
variables. These operating variables include powers consumed, power flows, 
power produced by generating units, load shed by demands, and voltage an¬ 
gles at buses (see reference [15] for a detail description of the formulation), 
n represents the maximum budget for investment in transmission expansion. 
Finally, f2(a:, d) dehnes the feasibility region for the operating variables y, as 
a function of investment decisions x and given realizations of the uncertain 
parameters d, as follows: 
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where A, B, F, G and K are matrices of constant parameters depen¬ 
dent on the network conhguration and element characteristics, Jeq selects 
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the components of y that are equal to the uncertain parameters (demands), 
and /ineq selects the components of y that are limited by the uncertain pa¬ 
rameters (i.e. maximum power generation and maximum load shedding). 
The first set of equality constraints correspond to constraints enforcing the 
power balance at every bus, the power flow through each line, and fixing the 
voltage angle of the reference bus. The first set of inequality constraints are 
associated with line flow limits, and limits on the voltage angles at every bus. 
Note that A, /x, a and are the dual variable vectors associated with each 
set of constraints, respectively. 

For a detailed physical interpretation of the mathematical formulation 
(l)-(3) we recommend the paper by [14]. 

2.1. Bi-level approach proposed by [15] 

[15], analogously to [12], proposed dealing with problem (l)-(3) by decom¬ 
posing and iteratively solving a subproblem and a master problem. The mas¬ 
ter variables correspond to x, i.e. the vector of first stage binary variables. 
Thus, for given values of these master variables, the subproblem corresponds 
to: 

Maximum Minimum y . (5) 

d E D y E Q{x, d) 

Considering the dual problem associated with the third-level, (5) results 
in the following single level maximization problem; 

Maximize yduai _ — (K — Fx^y, -I- (f{cx — ip) (6) 

d, X,y,a,(p 

subject to 


B^X-G^y + I^^a-lL^cp = b (7) 

y > 0 ( 8 ) 

‘F > 0, (9) 

d E D. (10) 


Regarding the uncertainty set D, [15] proposed the use of the following 
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set: 
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where df is the uncertain generation limit related to the ith. variable within 
vector d, df is the corresponding nominal value, df is the maximum positive 
distance from the nominal value that can take the random parameter, and 

is the set of indices of the generating units. Analogously, df, df, df 
and correspond to the same values but for demands. The only limitation 
introduced by this simplihed set is that uncertainty budgets and 
must be integer values, which in our opinion does not dwarf the benehts of 
robust optimization from the practical perspective. Alternatively, instead of 
uncertainty budgets associated with generation, demand and regions and 
r^, a unique uncertainty for each region T^, or for the system T could be 
used instead. 

Note that this set takes advantage of the simplihcation proposed by [17] 
because the uncertain parameters might only take either the nominal, or 
upper or lower limits of their uncertainty range, and the fact noticed by [14] 
that the worst realization of “nature” try to make generation capacities to 
be as lower as possible and demand loads as higher as possible. This set 
contains a reduced number of binary variables with respect to [12], which 
reduce the number of binary variables by half, increasing the percentage of 
reduction with respect to [13] and [14]. 

Subproblem (6)-(10) is a bilinear mathematical programming problem, 
which can be linearized and transformed into a mixed-integer linear mathe¬ 
matical programming problem at the expense of introducing binary variables 
associated with the uncertainty set (for more details, see the transformation 
of d'^{cy. — cp) introduced by [15]). 

The optimal solution of subproblem (6)-(10) provides the uncertain pa¬ 
rameter values (i(i) to construct primal cuts for the master problem, which 
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corresponds to the following optimization problem at iteration k\ 

Minimize c^x + 'j (17) 

x,y(iy, Wi = l,...,k-1 
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Note that the master problem, besides variable 7 , includes one operating 
variable vector for each realization of the uncertain parameters obtained 
from subproblem ( 6 )-( 10 ). 


3. Proposed solution technique 

This section justihes and explain in detail the proposed solution tech¬ 
nique. However, before focusing on the mathematical description of the so¬ 
lution method, we justify the type of uncertainty set used in this paper. 

3.1. Ellipsoidal uncertainty set 

Analogously to [3] and [ 6 ], we propose the following uncertainty set. Let 
us assume that uncertain parameter vector d = (d*^, has nominal or 

expected values d and variance-covariance matrix S, respectively. According 
to [3], the ellipsoidal uncertainty set can be written using the Mahalanobis 
distance as follows: 

£>(f)) = {(d-d)’'(E)-‘(d-d) <,?=}, (26) 

which constitutes the uncertainty set in (10). In contrast with the method 
proposed by [15], one unique parameter (3 controls the size and protection 
level of the ellipsoidal set. 


Note that only first and second order moments of the random parameters 
without lower and upper bounds are initially considered in this paper for the 
ellipsoidal uncertainty set. There are two reasons, firstly, the simplification 
of the process as it will be shown in next subsections, and secondly, because 
it is straightforward to include also the interval limitations (^d — d,d + ^ 
of the random parameters within the proposed solution method. 

In order to simplify the problem, the uncertainty set (26) is transformed 
using an affine mapping into a ball of radius /3, resulting in the following set 
of alternative constraints: 


d = d + Lz; \\z\\ < (3, (27) 

where L is the mapping matrix which can be obtained from Cholesky decom¬ 
position of variance-covariance matrix S = LL^, z represents a perturbation 
vector and || • || stands for Euclidean norm. Note that using also the fact no¬ 
ticed by [14] that the worst realization of “nature” try to make generation 
capacities to be as lower as possible and demand loads as higher as possible, 
we could reduce the uncertainty set by half including the following additional 
constraints: 


d^ < d^ (28) 

(29) 

which force maximum power productions and demand loads to take, respec¬ 
tively, values below and above their corresponding nominal values. 

3.2. Relationship of subproblem with First Order Second Moment methods 
Structural reliability techniques have been applied successfully to solve 
certain class of stochastic programming problems [18], and also robust opti¬ 
mization design problems [6]. Trying to establish an analogy between sub¬ 
problem (5) and the problem treated in [18], which attempts to optimize a 
given quantile of the random objective function, we could express the former 
as follows: 

/®(<i) = Minimum (b^ y), (30) 

y G 0.(x, d) 

however, no information about the worst realization of nature is included 
yet, for that reason it is necessary to use the optimal value df^^ of the un¬ 
certain parameters within the uncertainty set (27)-(29) which provide the 
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worst-minimum operational costs. It is obtained from solving the following 
problem; 

qp = maximum [/®(<i)] . (31) 

d G D{I3) 

Let assume that is known, in that case it is straightforward to cal¬ 
culate the optimal solution of subproblem (5) and values of the operating 
variables y* as follows: 

qp = Minimum y, (32) 

y G VL{x^ d) 


subject to 

d = d^^ (33) 

where 77 is the dual variable vector associated with constraint (32), which pro¬ 
vides information about how much the optimal objective function changes 
when uncertain parameters change around For those values of the opti¬ 
mal operating variables y*, we could try to hnd a new value of the uncertain 
parameter vector in the vicinity of d^'" but providing an even worst operating 
cost as follows: 


Maximize 

d 


subject to 


9/3 


dd 




T 


■ {d- 

- d^L) = qp + ri^{d - d^L) 

(34) 

d^ 

< 

d^ 

(36) 

d° 

> 

d , 

(36) 

d 

= 

d - 1 - Lz, 

(37) 

( 1 / 2 ) 

< 


(38) 


Since the ellipsoidal uncertainty set dehnes a convex region, the optimal 
solution d^'" must be located at the boundary of the ellipsoidal uncertainty 
set, i.e. constraint (38) must hold as an equality constraint. The optimal 
solution of the random parameters from problem (34)-(38) corresponds to 
d'^^, otherwise it would mean that it is possible to obtain another value of 
the random parameters inside the ellipsoidal set providing a worst value of the 
operational cost, so that d^'" would not correspond to the optimal solution, 
contradicting our initial hypothesis. This result might seem rather obvious. 
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however, this formulation (32)-(38) allows the proposal of a decomposition 
technique to look for the optimal solution of subproblem (5). 

Before describing the algorithm, let transform the formulation (34)-(38) 
into the following equivalent problem at the optimum (see reference [18] for 
the formal proof): 

P = Minimum 

d 

subject to 

d = 

qp + 

Note that the last constraint in (39) is the linear approximation of the 
equation: 

g{x,y*,d) = f{d) > q/s (40) 

which avoids using operating variables and solving (30) to get (40). This 
problem structure is analogous to that given by [18], however in that work 
the objective function depends explicitly on the uncertain parameters and 
there is no need to use a hrst order approximation. 

Problem (39) corresponds to the formulation associated with the struc¬ 
tural reliability method FOSM [see 19, 20, 21, 22, 23, 24, 25] used to cal¬ 
culate the probability of obtaining an operational cost higher than i.e. 
Prob(/®(£Z) > Qp). The random variables involved d belong to an u-dimensional 
space, which for given values of the decision x and operating y* variables 
can be divided into two regions with respect to the limit-state equation 
g{x,y*,d) = f{d) — qp, the undesirable {U) and desirable (V) regions. That 
is, 

U = {djg(x,y*,d) = f^(d) -q^ > 0 }, , . 

V = {dHg(x,y*,d) = f(d)-q 0 <0}. 

Note that once the qp associated with the worst operational costs within 
the uncertainty set is known, obtaining a greater value of that cost associated 
with another realization of the uncertain parameters outside the uncertainty 
set is an undesirable situation. 

In this context, the dependent normally distributed random variables d 
are transformed into independent standard normal random variables 2 that 
can be integrated straightforwardly using the linear Hasofer and Lind trans¬ 
formation [20], L is the Cholesky factorization of the variance-covariance 




d + Lz, 


dd 


{d - d“^) > q 0 . 


(39) 
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matrix, i.e. S = LL^. This factorization is always possible for positive def¬ 
inite matrices, and variance-covariance matrices must be positive definite by 
definition because their eigenvalues represent variances which must be always 
positive. This transformation is equivalent to the mapping transformation in 
(27). The optimal solution z* corresponds to the closest point to the origin 
located on the limit-state equation in the standard and independent normal 
random space, is the minimum distance so-called reliability index in the 
structural reliability scientific community, and d* = is the point of max¬ 
imum likelihood^ i.e. the actual values of the uncertain parameters located on 
the limit-state equation where the probability is higher, and it represents the 
most likely values of the random parameters that produce a non desirable 
operational cost. 

The hnal probability of obtaining an undesirable operational cost is re¬ 
lated to the reliability index by the relation (see [6] or [18] for more details): 

Prob(f(d)>g^) =$(-/?), (42) 

and the probability of a desirable operational cost is 

Prob {f{d) <qy)=l- $(-/?) = $(/3), (43) 

where 4>(-) is the cumulative distribution function of the standard normal 
random variables. This method provides the exact probability if the limit- 
state equation is linear in the standard normal random space, i.e., if the 
resulting limit-state distribution is normally distributed, which is the case if 
uncertain parameters are normally distributed. Nevertheless, it is an approx¬ 
imation in case parameters do not follow a multivariate normal distribution, 
which is usually the case becasue we only use expectations d and variance- 
covariance matrix S. According to expressions (42) and (43), qy corresponds 
to a upper quantile of the uncertain minimum operational cost random vari¬ 
able, or an approximation of the quantile, depending on whether the random 
parameters involved are normally distributed or not. 

From this perspective, subproblem (5) can be interpreted as finding the 
quantile of the uncertain minimum operating costs given by expression (43). 
According to [18] the solution method proposed in this paper attempts to 
minimize the Value-at-Risk of total transmission costs. 

Figure 1 shows a graphical illustration of the interpretation of subproblem 
5 under FOSM method. Panel (a) shows the circumferences associated with 
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different transformed joint probability density function contours for the bi- 
dimensional case (two random variables involved), including that associated 
with the ellipsoidal uncertainty set (gray dashed line) whose distance from the 
origin is equal to /9. It is also shown the operating cost contour equal to the 
quantile qp, which is tangent with respect the ellipsoidal set circumference at 
the point of maximum likelihood. Note that combination of random variables 
whose operating costs are higher than qg are in the undesirable region. If 
random variables are normally distributed, those operating cost contours 
are linear in the transformed space and the multidimensional problem can 
be reduced to one single dimension related to the standard normal random 
variable (see panel (b)). The probability of the operating costs to be inside 
the undesirable region, which requires a multidimensional integration can be 
calculated as $(—/?) and it corresponds to the gray shadow area in panel (b). 
This probability can also be represented in the space of the operating cost 
random variable (see panel (b)), where $(—/?) is equal to the shadow area 
associated with operating costs exceeding 



Figure 1: 3-D graphical interpretation of the FOSM method. 
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3.3. Proposed decomposition method 

Once the uncertainty set is properly defined and the relationship between 
subproblem (5) and FOSM methods has been established, we focus on our 
methodological proposal for solving problem (l)-(3). It is based on iteratively 
solving the following problems: 

Subproblem: Consisting on formulation (32)-(38), which is the equivalent 
with respect to (5) including the uncertainty set definition (27)-(29). 

Master problem: It consist of the minimization of the objective function 
(17) subject to constraints (18)-(25). 

This iterative scheme is known as “outer” because it tries to solve the 
complete TNEP problem. This outer iterative scheme have proved to be 
successful for solving TNEP problems using cardinality constrained uncer¬ 
tainty sets, however, the key for the success of the proposed methodology is in 
the resolution of the min-max subproblem. We propose a block coordinate 
descent algorithm [26] that takes advantage of sub-subproblems (32)-(33) 
and (34)-(38) to solve each of them optimizing the objective function over 
a subset of variables while the remaining variables are fixed to given val¬ 
ues. This process is iteratively solved until no further improvement of the 
objective function is achieved. This iterative scheme is known as “inner”. 
Convergence of this algorithm is only guaranteed for convex problems where 
each optimization sub-subproblem has a unique optimizer (strict convexity 
satisfies this requirement), which is the case for both sub-subproblems in this 
specific application. 

Before proceeding to explain how the subproblem can be solved, the pro¬ 
posed outer iterative scheme is described step by step in the following algo¬ 
rithm, which is analogous to that proposed by [14] and in [15]: 

Algorithm 3.1. (Outer iterative method). 

Input: Selection of uncertainty factor fd and the tolerance of the process e. 
These data are selected by the decision maker. 

Step 1: Initialization. Initialize the iteration counter to z/ = 1, and upper 
and lower bounds of the objective function = oo and = —oo. 
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Step 2: Solving the master problem at iteration u. Solve the master 
problem (17) subject to constraints (18)-(25). The result provides val¬ 
ues of the decision variables and 7 (y). Update the optimal objec¬ 
tive function lower bound Note that at the hrst 

iteration the optimal solution corresponds to the no investment case, 
alternatively, we could start with any other vector of decision variables. 

Step 3: Solving snbproblem at iteration v. For given values of the de¬ 
cision variables we calculate the worst operating costs within the 
uncertainty set or quantile obtaining also the corresponding un¬ 
certain parameters This is achieved by solving sub-subproblems 

(32)-(33) and (34)-(38). Update the optimal objective function upper 
bound = c^X(^y) + 

Step 4: Convergence checking. If < e go to Step 5, 

else update the iteration counter iz ^ u -\-l and continue in Step 2. 

Step 5: Ontpnt. The solution for a given tolerance corresponds to x* = 
X{u)- 


Next step is the description of the method to solve the subproblem, as 
pointed out previously it consists on a block coordinate descent algorithm 
that takes advantage of formulations (32)-(33) and (34)-(38). The key is to 
chose an initial value of the uncertain parameter vector at iteration k = 1, as 
candidate for maximum likelihood point d^. We recommend the selection of 
low values of generation capacities and high values of demand loads because 
we know that the optimal maximum likelihood point is around this area of 
the uncertainty set. Once this selection is done, sub-subproblem (32)-(33) is 
solved replacing constraint (33) by: 

d=d^^:r]^^y (44) 

The solution provides the optimal values of operating variables 
for that particular realization of the random parameters and the derivatives 
of the corresponding worst operational cost with respect to the uncertain 
parameters, i.e. '?7 (k)- Next, we try update the point of maximum likelihood 
in a attempt to hnd either a point inside the ellipsoidal uncertainty set, this 
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case occurs if the initial point is outside the ellipsoidal set, or a point in 
the boundary of the ellipsoidal set with an even worst operational cost. This 
is done solving sub-subproblem (34)-(38) replacing the objective function 
(34) by: 

Maximize + 77 J)) (d - djf . (45) 

d 

The advantage of using this sub-subproblem is that [6] provides the ana¬ 
lytical solution, which for this particular case is as follows: 

= d + (46) 

This inner” iterative process continues until no further improvement in 
the worst operational cost is achieved, and/or the maximum likelihood point 
does not change between consecutive iterations. The corresponding algo¬ 
rithm is summarized as follows: 

Algorithm 3.2. (Inner iterative method). 

Input: The same input arguments than the outer algorithm are used and 
also the actual values of the decision variables from the outer problem 
which are constants during this iterative process. 

Step 1: Initialization. Initialize the iteration counter to k = 1, and the 
initial values of the random parameters d^. 

Step 2: Determining the minimum operational cost at iteration k . 

For given values of the random parameters dJ^, obtain the operating 
variables providing the lowest operational costs for the actual iteration 
by solving problem (32)-(33) but replacing constraint (33) by (44). The 
result provides the optimal values of operating variables for 
that particular realization of the random parameters and the deriva¬ 
tives of the corresponding worst operational cost with respect to the 
uncertain parameters, i.e. '»7 (k)- 

Step 3: Convergence checking. If the iteration counter is equal to {k = 
1) go to Step 4i otherwise, proceed to check convergence, i.e., if — 
II < e go to Step 5, else continue in Step 4- 


16 




Step 4: Evaluating the new points of maximum likelihood. For given 
values of the operating variables we update the point of maximum 
likelihood using the analytical expression (46). Thus, updated repre¬ 
sentative values of the random parameters for the next iteration 
are obtained. The iteration counter must then be updated k ^ k + 1 
before proceeding to Step 2. 

Step 5: Output. The solution for a given tolerance corresponds to 

yU and 


3.4- Bounded random variables 

So far, the method developed in this paper deals with uncertainty coef¬ 
ficients by using the first and second probability moments associated with 
their probability distributions, and without assuming any particular type of 
distribution. This feature makes the method to be an approximation unless 
the random parameters truly follows a normal distribution. However, most 
of robust uncertainty sets proposed in the literature deal with random vari¬ 
ables within a given interval. This information could be easily incorporated 
in the proposed method by including the following constraint set into the 
subsubproblem (34)-(38): 

—d < d — d < d. (47) 

These constraints constitute bounded convex polytopes both in the origi¬ 
nal space associated with random model parameters d and also in the trans¬ 
formed space 2 . To solve this modified subproblem, the best strategy is as 
follows: 

1. To use the analytical solution (46) and check afterwards whether con¬ 
straints (47) hold, which means that the optimal solution has been 
achieved. 

2 . Otherwise, it means that the solution is not optimal according to those 
bounds, then the optimal solution is found by solving problem (34)- 
(38), including constraints (47) and removing the quadratic constraint 
(38). This results in a linear mathematical programming formulation. 

3. Finally, in case the solution from item 2 is also out of the ellipsoidal 
uncertainty set, the complete problem (34)-(38) including constraints 
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(47) must be solved. This is a linear mathematical programming prob¬ 
lem with one quadratic constraint that can be solve efficiently using 
interior point algorithms. 

Alternatively, after checking the analytical solution from item 1, the conic 
problem in item 3 could be used directly instead of the linear problem in item 

2 . 

4. Numerical case studies 

In this section, we present numerical experiments of our model proposal 
and its comparison with the method proposed by [15], which proved to be 
the most computationally efficient method to solve the TNEP problem using 
cardinality constraint uncertainty sets. We use the same examples provided 
by [15], i.e. the Garver system illustrative example [27], and two realistic 
case studies: the IEEE 24-bus system [28] and the IEEE 118-bus test system 
[29]. To avoid providing redundant information numerical data will not be 
displayed in this manuscript, those readers interested in that information can 
hnd it in [15] and references therein. 

All examples have been implemented and solved using GAMS [30, 31] 
and GPLEX 12, respectively, on a PG with four processors clocking at 2.39 
GHz and 3.2 GB of RAM memory. Note that computing times reported 
correspond to the sum of running times used by the solvers in order to solve 
the masters and subproblems until the hnal solution is attained. We use the 
GAMS command model. resusd. We use the same tolerance for all problems 
and equal to e = 10“®. 

4-1. Illustrative example. Garver system 

The proposed model is illustrated with the Garver 6-bus system, depicted 
in Figure 2. This system comprises 6 buses, 3 generators, 5 inelastic demands 
and 6 lines. Nominal values of generation capacities and demands and their 
offering and bidding prices can be found in [15]. The load-shedding cost is 
equal to the bid price of each demand. It is considered that a maximum of 
three lines can be installed between each pair of buses. Line data are obtained 
from Table I of [32], including construction costs, and the maximum available 
investment budget is €40 million. 

The investment return period of each line is considered to be 25 years, and 
the discount rate is 10%, resulting in an annual amortization rate of 10%. 
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-Existing lines 

-New lines null correlation 


New lines negative correlation 



Figure 2: Carver’s 6-bus test system. 


Since investment cost is annualized, the weighted factor a is equal to the 
number of hours of a year, i.e. 8760, obtaining an annualized load-shedding 
and power generation cost. 

Regarding the uncertainty set, the standard deviations of power gener¬ 
ation capacity and demand load are equal to 50% and 20% of their nomi¬ 
nal values, respectively, divided by 2.3263. Note that we have selected this 
standard deviation so that in case the random parameters are normally dis¬ 
tributed, the probability of those random parameters to be within the range 
dehned by [15] is equal to 0.99. Initially the random parameters are assumed 
to be independent. 

We have solved the robust TNEP problem for different reliability indexes 
(3. These values have been selected in an attempt to hnd similar objective 
functions with respect to those obtained in [15] for different combinations of 
the uncertainty budgets and T^. In order to obtain statistically sound 
conclusions about computing times, we have solved the problem 100 times 
and obtained mean and standard deviation executions times. Results from 
this numerical experiments are given in Table 1. 
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Table 1: Computational results for Carver’s 6-bus example 



[15] method 

Proposed method 


Optimal 
sol. (M€) 

H iter. 

Mean (s) 

Std. (s) 

P 

Optimal 
sol. (M€) 

H iter. 

(s) 

at (s) 

= 3 
pD = 5 

35505.31 

3 

0.644 

0.054 

6.58 

35505.31 

2 

0.16 

0.02 

D Q 

II II 

CO to 

25832.22 

4 

0.936 

0.094 

4.3 

25841.82 

3 

0.95 

0.05 

pG = 1 

r° = 2 

5861.92 

4 

1.011 

0.072 

2.6 

5315.16 

3 

0.44 

0.06 

D Q 

II II 

o o 

440.07 

2 

0.428 

0.065 

0 

440.07 

2 

0.17 

0.06 


From this table the following observations are pertinent: 

1. The proposed method is comparable in average computational time 
with respect to the method proposed by [15] for all cases, from the hrst 
and worst possible case corresponding to Soyster’s solution [1] up to the 
deterministic case when uncertainty budgets are null. Regarding the 
standard deviation of computing times the values are also comparable. 

2. Both methods require an small number of outer iterations. 

3. The worst possible solution associated with jS = 6.58 is the same than 
the analogous proposed by [15]. The reason is that all random parame¬ 
ters take their minimum and maximum possible values associated with 

2 ;j-values, i.e. equal to 2.3263. Note that 2.3263^ = 6.58. 

4. The deterministic solution is also the same because it corresponds to 
the case where all uncertain parameters are equal to their nominal 
values. 

5. For intermediate situations similar solutions have been found, which 
according to the reliability indexes (3 = 4.3 and {3 = 2.6 used and as¬ 
suming normally distributed parameters, they correspond to quantiles 
4>(/3) of 0.99999146 and 0.99533881. 

If we focus on the robust solution associated with the reliability index 
(3 = 4.3, the lines constructed are shown in Figure 2 in broken lines. At 
the point of maximum likelihood, only generation capacities of generators 
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at nodes 3 and 6, i.e those with higher capacities, reach their minimum 
values. The expansion investment is equal to 21.239M€. However, let assume 
that those power generators correspond to wind farm helds and their power 
productions are negatively correlated with coefficient —0.8, this makes the 
possibility of both generators at nodes 3 and 6 to be at their minimum 
capacities simultaneously remote. If we solve the robust problem under this 
new circumstance, the annualized optimal solution is reduced to 7739.721M€ 
with the same level of uncertainty, however, the cost of investment is higher 
38.616M€ because besides the additional lines constructed in the previous 
case, three additional lines joining buses 4 and 6 are constructed (shown in 
gray continuous lines in Figure 2). In this case, at the point of maximum 
likelihood, only generation capacities of generators at nodes 1 and 3 reach 
their minimum values. This result proves that the correlation structure of 
generation capacities associated with renewable and intermittent sources is 
crucial for the appropriate expansion of the network. 

Finally, we perform an additional simulation experiment using a reliabil¬ 
ity index [5 = $“^(0.9) = 1.28155 to check whether the optimal solution of 
the subproblem corresponds to the 0.9-quantile. Note that we have selected 
this level of protection because the maximum likelihood point at the optimal 
solution does not contain any limit of the random variables, in addition we 
divide the standard deviations of all random variables by 10 to make sure 
that we avoid sampling negative values of generations capacities and demand 
loads. This is important if we are interested in checking whether or not the 
optimal = 440.227M€ for this example corresponds to the 0.9-quantile. 
We also consider the negative correlation between generation capacities at 
nodes 3 and 6. The procedure is simple, we sample random values of gener¬ 
ation capacities and demand loads using the following expression: 

= d + (48) 

where are randomly generated samples from an independent standard 
normal distribution. For each realization of the uncertain parameters, prob¬ 
lem (32)-(33) replacing (33) by the expression d = is solved, which 
allows calculating the minimum operating costs g™ for those uncertain pa¬ 
rameters. Repeating the process 1000 times, the probability of not exceeding 
the optimal value = 440.227 is equal to the number of times the cost is 
lower or equal than g| divided by the total number of samples. The prob¬ 
ability obtained from Monte Carlo simulation is 0.907 which conhrms that 
the optimal cost corresponds to the 0.9-quantile. 
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Figure 3 shows the histogram associated with operational cost simula¬ 
tion results, the normal distribution htted and also the optimized quantile 
= 440.227. Note that the cost distribution is normal, because uncertain 
parameters were supposed to be normally distributed. For this reason the 
quantile minimization is exact. 

Garver simulation experiment 

O."" 

0.18 
0.16 
0.14 
0.12 

■g 0.1 

o 

Q 

0.08 
0.06 
0.04 
0.02 

Operational eost (Memo) 

Figure 3: Simulation experiment for Carver’s 6-bus test system. 


Garver simulation experiment 
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Althongh this resnlt is satisfactory, let us remind readers that this ex¬ 
act probabilities have been obtained nnder very specihc circumstances which 
barely happen in practice, for the rest of cases the solution is an approxima¬ 
tion of the quantile. 

4-2. IEEE 24-bus Reliability Test System 

The following case study described is based on the IEEE 24-bus Reliability 
Test System (RTS) [33], depicted in Fig. 4. The system comprises 24 buses, 
34 existing corridors which can accept a maximum of three equal lines and 
7 new corridors, 10 generating units and 17 loads. Data for lines in existing 
corridors are taken from [33], while line data for new corridors are obtained 
from Table I in [34]. Investment costs are €20 million, and analogously to 
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the Garver illustrative example, the investment return period of each line 
is considered to be 25 years, and the discount rate is 10%, resulting in an 
annual amortization rate of 10%. Data related to location of generators 



Figure 4: IEEE 24-bus reliability test system (RTS) 


and demands in the network, maximum power offered, prices for generating 
units, power bid and the bid prices for demands can be found in [15]. The 
load-shedding cost is equal to 100 times the bid price of each demand. 

Regarding the uncertainty set, power generation capacity can increase 
or decrease a maximum of 50% of their nominal values, while demand are 
allowed to change a maximum of 20%. Regarding the uncertainty set and 
analogously to the previous case, the standard deviations of power generation 
capacity and demand load are equal to 50% and 20% of their nominal values, 
respectively, divided by 2.3263. Note that we have selected this standard 
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deviation so that in case the random parameters are normally distributed, 
the probability of those random parameters to be within the range defined 
by [15] is equal to 0.99. Initially the random parameters are assumed to be 
independent. 

We have solved the robust TNEP problem for different reliability indexes 
l3. These values have been selected in an attempt to find similar objective 
functions with respect to those obtained in [15] for different combinations of 
the uncertainty budgets T*^ and T^. In order to obtain statistically sound 
conclusions about computing times, we have solved the problem 100 times 
and obtained mean and standard deviation executions times. Results from 
this numerical experiments are given in Table 2. 


Table 2: Computational results for IEEE 24 RTS case study without correlation 



[15] method 

Proposed method 


Optimal 
sol. (M€) 

(j iter. 

Mean (s) 

Std. (s) 

/3 

Optimal 
sol. (M€) 

tt iter. 

ht (s) 

crt (s) 

= 10 
pD = 17 

500752.14 

3 

1.51 

0.30 

12.09 

500752.14 

2 

1.34 

0.19 

= 7 
rD = 12 

454917.47 

3 

1.80 

0.39 

9 

455488.75 

3 

2.35 

0.71 

= 3 
r° = 5 

348941.40 

3 

1.12 

0.11 

4.5 

345282.84 

3 

0.78 

0.10 

= 0 
r° = 0 

219161.52 

2 

0.414 

0.08 

0 

219161.52 

2 

0.28 

0.06 


From Table 2 we can extract the same conclusions than in the previous 
illustrative example: 

1. The proposed method is comparable in average computational time 
with respect to the method proposed by [15] for all cases, from the first 
and worst possible case corresponding to Soyster’s solution [1] up to the 
deterministic case when uncertainty budgets are null. Regarding the 
standard deviation of computing times the values are also comparable. 

2. Both methods require an small number if outer iterations. 

3. The worst possible solution associated with /3 = 12.09 is the same than 
the analogous proposed by [15]. The reason is that all random parame- 


24 

















ters take their minimum and maximum possible values associated with 
2 ;j-values, i.e. equal to 2.3263. Note that 2.3263^ = 12.09. 

4. The deterministic solution is also the same because it corresponds to 
the case where all uncertain parameters are equal to their nominal 
values. 

5. For intermediate situations similar solutions have been found. For 
instance, the optimal solution associated with the reliability index 
(3 = 4.5, and assuming normally distributed parameters, corresponds 
to quantile <l>(/3) of 0.999996602326875. 

If we focus on the robust solution associated with the reliability index 
(3 = 9, the lines constructed are shown in Figure 4 in broken lines. Only three 
new lines are selected, from bus 3 to 9, from 11 to 14 and from 14 to 16. The 
investment cost is 14.307M€. At the point of maximum likelihood, there 
are several generation capacities reaching their minimum possible values, for 
instance generators 3 and 4 at buses 7 and 13, respectively. If we assume 
that those power generators correspond to wind farm helds and their power 
productions are negatively correlated with coefficient —0.8, this makes the 
possibility of both generators at nodes 3 and 4 to be at their minimum capac¬ 
ities simultaneously remote. If we solve the robust problem under this new 
circumstance, the annualized optimal solution is reduced to 432351.604M€ 
with the same level of uncertainty, however, the cost of investment is higher 
16.606M€ and instead of constructing the lines from the previous case, only 
two lines from bus 3 to 24 and from 15 to 24 are selected as new lines 
(shown in gray continuous lines in Figure 4). This result also conhrms that 
the correlation structure of generation capacities associated with renewable 
and intermittent sources is important for the appropriate expansion of the 
network. 

Finally, we perform an additional simulation experiment using a reliability 
index (3 = $“^(0.9) = 1.28155 to check whether the optimal solution of the 
subproblem corresponds to the 0.9-quantile. Note that we have selected this 
level of protection because the maximum likelihood point at the optimal 
solution does not contain any limit of the random variables, which is the 
case with the standard deviation used in the example. 

This is important if we are interested in checking whether or not the op¬ 
timal q*^ = 251719.780M€ for this example corresponds to the 0.9-quantile. 
We also consider the negative correlation between generation capacities at 
nodes 7 and 13. The procedure is analogous to that used in the previous 
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illustrative example. Repeating the process 1000 times, the probability of 
not exceeding the optimal value = 251719.780 is equal to the number of 
times the cost is lower or equal than divided by the total number of sam¬ 
ples. The probability obtained from Monte Carlo simulation is 0.889 which 
conhrms that the optimal cost corresponds to the 0.9-quantile. 

Figure 5 shows the histogram associated with operational cost simulation 
results, the normal distribution htted and also the optimized quantile q*^ = 
251719.780. Note that the cost distribution is normal, because uncertain 
parameters were supposed to be normally distributed. For this reason the 
quantile maximization is exact. 


IEEE 24 RTS simulation experiment 



Figure 5: Simulation experiment for IEEE 24 RTS test system. 


Although this result is satisfactory, let us remind readers that this ex¬ 
act probabilities have been obtained under very specihc circumstances which 
barely happen in practice, for the rest of cases the solution is an approxima¬ 
tion of the quantile. 
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4-3. IEEE 118-bus test system 

Finally, we run additional computational tests using the IEEE 118-bus 
test system [29]. The system comprises 118 buses, 186 existing lines, 54 
generating units and 91 loads. In addition, it is possible to construct up to 
61 additional lines to duplicate each one of the following existing lines: 8, 
12, 23, 32, 38, 41, 51, 68, 78, 96, 104, 118, 119, 121, 125, 129, 134, 159, 7, 
9, 36, 117, 71, 131, 133, 147, 103, 65, 144, 168, 4, 13, 132, 69, 66, 67, 5, 
89, 29, 167, 145, 70, 42, 90, 16, 174, 98, 99, 185, 93, 94, 128, 164, 97, 153, 
146, 116, 163, 31, 92, 130. Data for lines in existing corridors are taken from 
[29]. Investment costs are €100 million, and analogously to both previous 
examples, the investment return period of each line is considered to be 25 
years, and the discount rate is 10%, resulting in an annual amortization rate 
of 10%. 

Data for generation capacities and demand loads are given in [15]. The 
load-shedding cost is equal to 1.2 times the bid price of each demand. 

Regarding the uncertainty set, the standard deviations of power gener¬ 
ation capacity and demand load are equal to 50% and 50% of their nomi¬ 
nal values, respectively, divided by 2.3263. Note that we have selected this 
standard deviation so that in case the random parameters are normally dis¬ 
tributed, the probability of those random parameters to be within the range 
dehned by [15] is equal to 0.99. Random parameters are assumed to be 
independent. 

We have solved the robust TNEP problem for different reliability indexes 
jS. These values have been selected in an attempt to hnd similar objective 
functions with respect to those obtained in [15] for different combinations of 
the uncertainty budgets T*^ and T^. In order to obtain statistically sound 
conclusions about computing times, we have solved the problem 100 times 
and obtained mean and standard deviation executions times. Results from 
this numerical experiments are given in Table 3. 

From Table 3 the same observations than in previous examples are ob¬ 
tained. However, in this case, the worst possible solution is associated with 

P = 28.02 because 452.3263^ = 28.02. Note that computational 

times are below 10 seconds, which prove the method to be applicable in real 
world energy transmission networks. 
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Table 3: Computational results for IEEE 118-bus test system 



[15] method 

Proposed method 


Optimal 
sol. (M€) 

tj iter. 

Mean (s) 

Std. (s) 

/3 

Optimal 
sol. (M€) 

tt iter. 

(s) 

(s) 

= 54 
r° = 91 

31994.36 

4 

7.73 

0.80 

28.02 

31994.36 

2 

0.42 

0.25 

= 35 

r° = 60 

30032.93 

5 

53.46 

60.36 

21.0 

29827.30 

4 

1.84 

0.49 

= 15 

= 20 

23352.97 

4 

24.26 

26.05 

9.5 

21828.82 

5 

4.60 

0.23 

= 0 
rD = 0 

13929.23 

2 

0.56 

0.07 

0 

13929.23 

3 

2.65 

0.34 


5. Conclusions 

This paper extends the applicability of method [6], based on first-order 
second-moment methods from structural reliability and decomposition tech¬ 
niques, into the adaptive robust optimization problem associated with trans¬ 
mission network expansion planning. The method is specially suitable for 
problems where first and second order moments of the probability distribu¬ 
tions of the uncertain parameters involved are available. 

The proposed method has the following advantages: 

1. The feasibility region defined by (27) is convex and only depends on 
variables d. 

2. The resulting sub-subproblem (32)-(33) for fixed values of the uncertain 
parameters is linear and it does not content binary variables, in con¬ 
trast with respect to TNEP problems based on cardinality constrained 
uncertainty sets. 

3. The sub-subproblem (34)-(38) resulting from fixing the operating vari¬ 
ables to given values has analytical solution. In case additional bounds 
on uncertain parameters are included, the worst possible scenario re¬ 
quires solving a linear programming problem with one quadratic con¬ 
straint, which can be solved efficiently using interior point algorithms. 

4. The outer iterative method between subproblem and master problem 
provides upper and lower bounds of the objective function, which al- 
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lows checking if the problem is being solved appropriately dnring the 
iterative process. Besides, it converges in a small number of iterations. 

5. This method allows including the correlation structure of the random 
variables involved easily, this is very important specially for wind power 
farms at different sites and demand loads, which might have significant 
correlations affecting the optimal expansion strategy. 

6. The interpretation based on structural reliability allows selecting the 
level of conservationism in terms of approximate probabilities of ex¬ 
ceeding the design operational costs, because this parameter allows 
selecting the quantile of the cost function to be optimized. This strat¬ 
egy is equivalent to the value-at-risk strategy widely used in financial 
and stochastic applications. 

In summary, this paper proposes a new decomposition algorithm to at¬ 
tain the exact solution of the TNEP problem derived from using a two-stage 
adaptive robust strategy with ellipsoidal uncertainty sets. In addition, this 
procedure set the basis for alternative and more sophisticated robust uncer¬ 
tainty sets based on reliability methods. This would allow incorporating not 
only expected values and the correlation structure of uncertain parameters 
but appropriate marginal distributions of those parameters. This is a subject 
for further research. 
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